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Abstract 

We study the phase behavior and the interfacial tension of the screened Coulomb (Yukawa) 
restricted primitive model (YRPM) of oppositely charged hard spheres with diameter a using 
Monte Carlo simulations. We determine the gas-liquid and gas-solid phase transition using free en- 
ergy calculations and grand-canonical Monte Carlo simulations for varying inverse Debye screening 
length k. We find that the gas-liquid phase separation is stable for kg < 4, and that the critical 
temperature decreases upon increasing the screening of the interaction (decreasing the range of 
the interaction). In addition, we determine the gas- liquid interfacial tension using grand-canonical 
Monte Carlo simulations. The interfacial tension decreases upon increasing the range of the inter- 
action. In particular, we find that simple scaling can be used to relate the interfacial tension of 
the YRPM to that of the restricted primitive model, where particles interact with bare Coulomb 
interactions. 
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I. INTRODUCTION 



Coulombic interactions are important in a wide variety of physical systems such as elec- 
trolytes, molten salts, plasmas, colloidal suspensions, micelles, microemulsions, and liquid 
metals. The screened Coulomb (Yukawa) potential arises naturally for charged particles in 
the presence of a screening distribution of microions. The phase behavior of a pure sys- 
tem of hard spheres interacting with screened Coulomb potentials has been well-studied 
and the phase diagram displays stable fluid, face-centered-cubic (fee), and body-centered- 
cubic (bec) crystal phases.— In this paper, we study a binary fluid of oppositely charged 
particles using computer simulations. While the phase diagram of the restricted primitive 
model (RPM), consisting of a binary mixture of equally sized hard spheres carrying oppo- 
site charges of equal magnitude, and interacting with bare Coulombic interactions, has been 
widely studied,— i 6 i 7 i 8 i 9 i 10 i n i 12 there is little information available on the phase diagram of the 
Yukawa restricted primitive model (YRPM), where the hard spheres of diameter a interact 
with screened Coulomb potentials Uij = ±ecexp[— «(ry — cr)]/r,j, with the distance be- 
tween particles i and j, e the contact value, and k the Debye screening parameter (inverse of 
the Debye screening length). Recently, Hynninen et al. 13 determined the full phase diagram 
of the YRPM for a screening parameter k<j = 6. At high temperatures, the system behaves 
like a pure hard-sphere system, with a transition between a fluid and a substitutionally 
disordered fee phase, where the opposite charges are distributed randomly on a fee lattice. 
At lower temperatures, a dilute gas phase coexists with a high density CsCl solid phase, and 
the gas-liquid transition is metastable with respect to freezing. At high densities, various 
solid-solid transitions appear, e.g., a transition from CsCl to CuAu and from CuAu to the 
tetragonal phase. Overall the system exhibits a phase behavior in striking similarity with 
the RPM phase diagram , 7 i 8 i 9 i 10 i 13 which displays a fluid-disordered fee transition at high 
temperatures, a stable gas-liquid transition at low temperatures, and a fluid-solid transition 
at higher densities. Since the RPM is the limit of the YRPM for na — > 0, we expect a 
crossover from a metastable to a stable gas-liquid transition for < kg < 6. 

The gas-liquid transition for a similar model with pair potential Uij = eexp[— tzfrij — a)], 
has been studied using computer simulations . 14 ' 1 ^ 16 This interaction potential differs by 
a factor a/r^ from our model. The factor a/r^ is of the order of unity, and we expect 
the effect on the phase behavior to be small. This allows us to compare the results of 
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Refs . 1 1 with the results of the YRPM. In particular, Caballero et al.— investigated 
the critical temperature as a function of the screening parameter. In a later paper,— the 
stability of the gas-liquid separation with respect to the gas-solid transition was estimated 
by computing the melting density of the CsCl structure. This technique overestimate the 
stability of the gas-liquid binodal, with respect to our free energy calculations. We will 
discuss the relationship between our results and those of Ref.— in more detail in Sec. IV. 

An experimental realization of the YRPM is provided by charge-stabilized colloidal sus- 
pensions. Recently, it was shown experimentally that the charge on the colloids can be 
tuned in such a way that oppositely charged colloids can form large equilibrium ionic col- 
loidal crystals . 1 ?i 1 Ti 1 ?i 1 ? Experiments, theory, and simulations based on screened Coulomb 
interactions are in good agreement. The system studied in Refs . 1 ^ 1 ? had a Debye screening 
parameter na ~ 7, and a gas-liquid phase separation was not observed. 

The critical temperature and structure of the YRPM has also been studied using integral 
equation theory— but no information on the stability of the gas-liquid transition with respect 
to freezing has been given. 

The purpose of this paper is to determine the value of the screening parameter kg at 
which the gas-liquid transition becomes stable for oppositely charged colloids. To this end, 
we perform Monte Carlo simulations to compute the Helmholtz free energies of the fluid 
and solid phases of the YRPM. We also study the dependence of the critical parameters 
and the gas-liquid interfacial tension on the interaction range. The critical parameters and 
the values of the interfacial tension are calculated using histogram reweighting methods and 
grand-canonical Monte Carlo simulations. The paper is organized as follows. In Sec. II, 
we describe the model, and we present the simulation methods in Sec. III. The results are 
discussed in Sec. IV and we end with some concluding remarks in Sec. V. 

II. MODEL 

We investigate the Yukawa restricted primitive model (YRPM) consisting of N spherical 
particles with a hard-core diameter a in a volume V. Half of the spheres carry a positive 
charge and the other half a negative charge of the same magnitude. The pair interaction 
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reads 

oo 

ey exp[— n(vij — a 



(3u(ri 



oo < a 



k B T rij/a 
otherwise 



a < < r cut , (1) 



where is the distance between spheres i and j, k the screening parameter, (3 = 1/ksT the 
inverse temperature with k B the Boltzmann constant and T the temperature, and ey the 
contact value of the potential. The cut-off value is r cut = 3.6a. The interaction is attractive 
for oppositely charged spheres, and repulsive for like-charged spheres. We define a reduced 
temperature T Y = ksT/ey and measure particle density in terms of the packing fraction 
v = (na 3 /6)N/V. 

According to the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory , 21 ' 22 the effective 
pair potential between two charged spheres carrying the same number Z of elementary 
charges e suspended in a sea of salt ions with density p s is given by Eq. (P) with a contact 
value 

k B T (l + Ka/2) 2 a 1 ' 

The Debye screening parameter reads k = y/EnXB/h, where As = e 2 /e s fcsT is the Bjerrum 
length and e s is the dielectric constant of the solvent. It must be noted that, more recent 
theories on same charged colloidal spheres suspended in a sea of salt ions yield potentials 
of the form of Eq. (P), but with screening parameters that depend on the charged colloid 
concentration . 23 ' 24 ! 25 ' 26 ! 2 i 28 i 29 However, the exact functional form is yet unknown 30 and dif- 
ferent theories predict varying functional forms. Furthermore, the DLVO theory was not 
originally derived for oppositely charged spheres, but it can be extended using the linear 
superposition approximation (LSA) to obtain the potential given by Eqs. ((TJ and (j2J)- 31 The 
extended DLVO theory has been shown to give good agreement with Poisson-Boltzmann 3 ^ 
and primitive model calculations 33 at small k<t, justifying the use potential in Eq. (JIJ) with 
the contact value given by Eq. (J2J). We will refer to the DLVO theory extended by the LSA 
simply as the DLVO theory. To facilitate the comparison between the results of the DLVO 
theory for different screening lengths k, we define a reduced temperature Tq = <j/Z 2 \b that 
is independent of k and equal to the the definition of the reduced temperature of the RPM. 
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III. SIMULATION METHODS 



In order to determine the stable phase for a given state point, we compute the Helmholtz 
free energy as a function of 77 and T y . As the free energy cannot be measured directly 
in a Monte Carlo simulation, we use thermodynamic integratio n 34 ' 3 ^ 3 ^ to relate the free 
energy of the YRPM system to that of a reference system, whose free energy is known. 
In the thermodynamic integration of the fluid phase, we use the hard-sphere fluid as the 
reference state, whereas in the solid phase, the reference state is the Einstein crystal. We use 
a 10-point Gaussian quadrature for the numerical integrations and the ensemble averages 
are calculated from runs with 40000 MC cycles (attempts to displace each particle once), 
after first equilibrating the system during 20000 MC cycles. Employing a common tangent 
construction on the fluid free energy density curves as a function of rj, we find the points of 
tangency that correspond to the densities of the coexisting gas and liquid phase. A similar 
common tangent construction is used to determine the coexistence between the fluid and 
solid phases, and to check whether the gas-liquid separation is stable with respect to the 
fluid-solid phase coexistence. In addition, we perform a more detailed study of the gas-liquid 
binodal using methods based on histogram reweighting. To this end, we employ grand- 
canonical Monte Carlo simulations with successive umbrella sampling 37 to overcome the 
free energy barrier between the liquid and gas phase. In the successive umbrella sampling 
method, the probability P(N)\ z+jX _ of having N particles at fugacity z = z + = z_ in a 
volume V = L 3 is obtained by sampling successively 'windows' of particle numbers. In 
each window, the number of spheres N is allowed to fluctuate by one particle, i.e., between 
and 1 in the first window, 1 and 2 in the second window, etc. We choose at random 
whether to make an attempt to insert or to remove a particle such that, on average, the 
system is charge neutral. The sampling of the probability ratio P(N)/P(N + 1) is done, in 
each window, until the difference between two successive samplings of the probability ratio 
is smaller than 10~ 3 . At phase coexistence, the (normalized) distribution function P(N)\ Z 
becomes bimodal, with two separate peaks of equal area for the liquid and gas phase. To 
determine phase coexistence we calculate the average number of particles 




(3) 
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Subsequently, we use the histogram reweighting technique^ 8 - to determine the fugacity z' for 
which the equal area rule 

(AT) poo 

P(N)\ z ,dN= / P(N)\ z ,dN, (4) 

J(N) 

which is the condition for phase coexistence, is satisfied. 

The gas-liquid interfacial tension 7i g for a finite system with volume V = L 3 is obtained 
from P(N)\ z i at coexistence: 



2L 2 



ln | gL) + p(KJ j _ ln(P(iVmin )) 



(5) 



where P(iV^ ax ) and P(A^ ax ) are the maxima of the gas and liquid peaks, respectively, and 
P{N m i n ) is the minimum between the two peaks. We determine the bulk interfacial tension 
7ig by performing simulations for a range of system sizes and by extrapolating the results to 
the infinite system size using the relatio n?^ 4 ^ 41 

_ xlnL InA 

where A and x are generally unknown. The finite size scaling is performed using the simu- 
lation results for box lengths L/cr =8, 10, 12, and 14. 

The critical temperature T cr and critical packing fraction determined by fitting 

the scaling law 

r )i- r l9 = h{Tcr - T) - 32 , (7) 
and the law of rectilinear diameters 

= Vcr + h(T cr - T), (8) 

to the simulation results for the gas (r] g ) and liquid (r/i) packing fractions, where fi and f2 
are fitting parameters. 



IV. RESULTS 



We compute the phase diagram using thermodynamic integration and grand-canonical 
Monte Carlo simulations for screening parameters k<j=2, 3, 4, 4.5, and 6. In Fig. d we 
show the resulting phase diagrams in the (rj, Ty) plane, together with the na=0 phase 
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FIG. 1: Phase diagrams of the YRPM in the reduced temperature T Y - packing fraction rj repre- 
sentation for varying Debye screening parameter (a) na = from Refs i 5 i 6 i 10 ; (b) ko = 2; (c) kct=3; 
(d) no = 4; (e) k<t = 4.5; (f) no = 6. The squares are the results from the free energy calculations 
and the circles are the results from the grand-canonical Monte Carlo simulations. F and S denote 
the stable fluid and solid (CsCl) phase. F + S and F + F denote, respectively, stable fluid-solid 
and (meta)stable gas-liquid coexistence region (the shaded regions are metastable). The dashed 
line in (f ) indicates the results of the model used by Caballero and Puertas M . The lines are a guide 
to the eye. Tie lines (not shown) are horizontal. 

diagram from Refs.-^ (fluid-solid) and 5 ' 6 (gas-liquid). The squares denote the results from 
the free energy calculations, and the circles represent the gas-liquid binodal obtained from 
the grand-canonical Monte Carlo simulations. We find good agreement between both results. 
The shaded areas in Figs, ^e)-(f) represent the metastable gas-liquid regions for screening 
parameters k<j = 6 and 4.5. For smaller screening parameters, Ka=4, 3, and 2, the gas-liquid 
transition is stable, and the phase diagram resembles that of a simple fluid. At sufficiently 
low temperatures, a gas-liquid phase separation (metastable for screening parameters kg = 6 
and 4.5) appears at low densities and a fluid-solid transition at high density. At the triple 
point, the gas, liquid, and the solid phase are in coexistence, while at the critical point, the 
gas and the liquid phase have the same density. At temperatures below the triple point, a 
dilute gas coexists with a high density solid, and at temperatures higher than the critical 
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temperature, a fluid coexists with a solid phase. Figure [T] shows that the region of stable 
liquid phase increases upon increasing the range of the interaction, i.e., decreasing kg. For 
simple fluids with short-range attractive Yukawa interactions, square-well attractions, and 
depletion attractions, the relationship between the range of the attractive interactions and 
the stability of the gas-liquid transition has been well-studied by computer simulations, 
density functional calculations, and integral equation theories . 42 i 43 i 44 i 45 i 46 i 47 i 48 These studies 
show that the minimum range of attractions required for a stable gas-liquid transition is 
about one sixth of the range of the repulsions. 

Our results show that the gas-liquid coexistence is stable for kg < 4, and therefore 
contradict the findings of Ref.- 1 ^, where an estimate kg < 25 was given. In this comparison, 
we have to keep in mind that the pair potential used in Ref-^ did not include the factor 1/r, 
which we include. To study the effect of this factor, we repeated our free energy calculations 
using the pair potential of Ref.— for screening length kg — 6. The results for this model 
are presented in Fig. WIS) with a dashed line. We find that there is no qualitative difference 
between the two models; both predict a metastable gas-liquid transition for screening length 
kg = 6. In Ref.^&, the stability of the gas-liquid transition was determined from the cross- 
over between the freezing line and the liquid branch of the gas-liquid binodal. The cross-over 
point was recognized as the gas-liquid-solid triple point. When the triple point was at a 
lower temperature than the critical point, the gas-liquid phase separation was considered 
stable. Since the CsCl structure melts in the middle of the broad gas-solid coexistence, the 
criteria used in Ref.^ would indicate a stable gas-liquid separation, whereas our free energy 
calculations show that it is, in fact, metastable. We argue that the computation of the 
melting line cannot be used to determine the stability of gas-liquid transition with respect 
to a broad gas-solid coexistence region. 

It is interesting to note that previous simulation studies of a one-component hard-core at- 
tractive Yukawa fluid predict a stable gas-liquid transition for kg = 3.9, while it is metastable 
for kg = 7 , 42 ^ which compares well with our results. 

Table I summarizes the critical temperatures T Y cr and critical packing fractions r] cv as 
found from our simulations and from the generalized mean spherical approximation (GMSA) 
theory^ for different values of kg. 

In Figure El we analyze the effect of the cutoff value of equation (JTJ) on the liquid-gas 
binodal. Huge deviations are expected as the interaction becomes longer ranged. We used 
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FIG. 2: Binodals of the YRPM for screening parameter kg=2, and cutoff values r cut /cx=2.3, 3.6, 
and 5 in the (rj, Tt) representation. Statistical errors are of the order of the symbol size. The lines 
are a guide to the eye. 
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TABLE I: Critical temperatures T Y cr and packing fractions r] CT for the YRPM for different values 
of the Debye screening parameter kg. The GMSA theory data is from Rcf. 20 . 

three different cutoff values r cut /<7 = 2.3,3.6,5 for the calculation of the liquid-gas binodal 
for our longest range interaction (kg = 2). The binodals for cutoff values r cut /cr = 3.6, and 5 
are equivalent within the statistical accuracy, thereby justifying the choice of a cutoff value 
r cut/0" = 3.6 in all subsequent calculations. 



9 



0.2 
0.15 
** 0.1 

< 

0.05, 

°0 0.1 0.2 0.3 0.4 

FIG. 3: Binodals and the critical points of the YRPM for screening parameters na=0, 2, 3, and 
4, in the (rj,Ty) representation. Crosses denote the location of the critical points. The binodal for 
k<7=0 (RPM) is taken from Ref£ and the critical point from Ref4^. The lines are a guide to the 
eye. 

Figure El shows the (stable) gas-liquid binodals for na=2, 3, and 4, in the (T),Ty) repre- 
sentation. Table I and Fig. 01 show that, for < kg < 6, the reduced critical temperature 
T Y cr and the critical packing fraction r] CT decrease for increasing range of the interaction, 
i.e., decreasing kg, in agreement with the findings of Caballero et al.— . The non-monotonic 
behavior of the critical temperature as a function of kg that was reported in Ref.-^ for screen- 
ing parameters kg > 10 is in the region where we claim the gas-liquid phase separation to 
be metastable. As can be seen from Table I, the GMSA theory predicts a non-monotonic 
behavior of the critical temperature as a function of kg. Comparing the theoretical results 
with our simulations, we observe that the GMSA theory overestimates the critical temper- 
ature for kg < 6 and underestimates it at kg = 6. On the other hand, the GMSA theory 
underestimates the critical packing fraction r] cr at kg = 0, but agrees reasonably well with 
our simulation results for kg > 3. 

Figure Ufa) shows the binodals of the YRPM and the RPM in the corresponding state 
representation, where the reduced temperature T Y is scaled with the critical temperature 
T Y cr , and the packing fraction r\ is scaled with the critical packing fraction i] cr . We see that 
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FIG. 4: (a) Binodals of the YRPM and the RPM in the corresponding state representation for 
Debye screening parameters kg = 2, 3, and 4. The reduced temperature T Y is scaled with the 
critical temperature Ty cr and the packing fraction r/ is scaled with the critical packing fraction rj CT . 
The binodal for kg=0 (RPM) is from Ref£ and the critical point from Ref.4£. (b) The difference 
in coexisting packing fractions Ar/i g = rji — r\ g is plotted against 1 — T y /T Y cr . 



the binodals do not collapse on a single master-curve, but instead, the RPM binodal (where 
kg = 0) differs considerably from the YRPM binodals (where kg = 2, 3, and 4). This finding 
is in agreement with the prediction of the GMSA theory— In Fig. Efb), we plot the width 
of the gas-liquid separation, Ai]\ g = r/i — r] 9 , as a function of 1 — T Y /T YjCr . We see that for 
a fixed 1 — T Y /T Ycr , the width of the gas-liquid separation decreases with increasing range 
of the interaction, resulting in a smaller density gap between the coexisting liquid and gas 
phase for longer-ranged interactions. 

Figure 03 shows the gas-liquid interfacial tension, scaled with the contact value energy 
ey, for different values of the screening parameter kg. For comparison, we also show the 
interfacial tension of the RPM from Ref.-^. As can be seen from Fig. the value of the 
dimensionless interfacial tension increases with increasing kg. This can be understood on the 
basis of Fig. Efb), which shows that, with increasing kg, the density gap of the coexistence 
region increases, meaning that the interfacial tension increases. The inset of figure El shows a 
log-log plot of 7i g <7 2 /e Y versus 1 — T Y /T Ycr in the vicinity of the critical point, which can be 
used to extract an estimate for the critical exponent of the correlation length v. We found 
1v ~ 1.1 for all screening parameters kg, by performing a linear fit on the data, which differs 
from the Ising model result 2v = 1.32,^2 and from the accepted value of 2u = 1.26. 51 ' 52 The 
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FIG. 5: Dimensionless gas-liquid interfacial tension 7i g cr 2 /eY as a function of 1 — T Y CT /T Y - The 
data for kg = is from Refi^. Inset: log- log plot of the dimensionless interfacial tension in the 
vicinity of the critical point. The dashed line is the theoretical prediction 51 ' 52 7 ~ (T Y — TZ^ 1 ~" 



value of the correlation length is very sensitive to the extrapolated surface tension to infinite 
system sizes. In order to improve the statistical accuracy of the simulations, larger system 
sizes, as well as longer runs are needed, especially close to the critical point. Nevertheless, 
our results are compatible within the simulation error, with the theoretical prediction of the 
correlation length (dashed line in the inset of figure EJ). 

We now interpret our results in view of the DLVO theory. In the DLVO theory, the 
contact value ey, and hence the reduced temperature T Y = hsT/ey, depend on the salt 
concentration p s through the screening parameter k<j, see Eq. (J2J). In Fig. E(a), we plot the 
gas-liquid binodals and critical points for kct=0, 2, 3, and 4, using the reduced temperature 
Tq = <j/Z 2 \b that does not depend on ko. As can be seen from Fig. Efa), the reduced 
critical temperature Tq cr decreases with increasing k<j, or salt concentration. This means 
that, at a fixed Tq and at a statepoint inside the gas-liquid coexistence region, adding salt 
decreases the density difference between the gas and the liquid phases, until, at the critical 
salt concentration, the density difference disappears. This finding could be confirmed by 
performing simulations with explicit co- and counterions , 30 ' 53 and could be used to experi- 
mentally test the validity of the DLVO theory for oppositely charged colloids. Figure Efb) 
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FIG. 6: (a) Binodals of the YRPM, for Debye screening parameters na =0, 2, 3, and 4, in the 
reduced temperature Tq and the packing fraction rj representation. The binodal for kct=0 (RPM) 
is from Ref£ and the critical point from Ref. 49 . (b) Dimensionless gas-liquid interfacial tension 
7i g o- 2 /ec as a function of 1 - T*/T* fiT . Remember that T*/T* cr = T*/T*^ CT . The data for kg = 
is from Ref. 50 . 



shows the interfacial tension scaled with the contact value ec = z2kB J XB ; and as can be 
seen, the interfacial tensions for different ko collapse to a single line. This suggests that 
the interfacial tension is determined solely by the contact value and not by the range of the 
interaction. 



V. CONCLUSIONS 

We have used a combination of MC free energy calculations and grand-canonical MC 
simulations to determine the stability and the interfacial tension of the gas-liquid phase sep- 
aration in a binary mixture of oppositely charged hard spheres, which interact via screened- 
Coulomb (Yukawa) potentials. We find that the gas-liquid coexistence is stable with respect 
to gas-solid coexistence for values of the screening parameter kg < 4. This value is similar 
to what is found for the single component attractive Yukawa model," 4 ^ where the gas-liquid 
transition is stable at ko = 4 and metastable at ko = 7. 

We have studied the dependence of the critical temperature as a function of the range of 
the Yukawa interaction. If the contact value of the interaction potential does not depend 
on the screening length, it is possible to define a reduced critical temperature simply as the 
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inverse of the Yukawa contact value. With this definition, the reduced critical temperature 
decreases upon increasing the range of the interaction, which is in agreement with Ref.— . 

We have related the Yukawa restricted primitive model (YRPM) to the DLVO the- 
ory, which was recently used to explain experimental results on oppositely charged 
colloids . 1 ? 11 ^ 1 ^ 1 ^ 4 The DLVO theory predicts a contact value that depends on the screening 
length. Thus, in order to facilitate the comparison between the results for different screen- 
ing lengths, we define a temperature scale that is independent of the screening length. The 
natural choice is the reduced temperature of the RPM, which is the limit of zero screening 
length of the DLVO theory. With this definition, the reduced critical temperature decreases 
upon increasing the range of the interaction. This means that upon adding salt to a system 
at fixed temperature and at a statepoint in the gas-liquid coexistence region, the density 
difference between the gas and liquid phases decreases, and finally disappears at the critical 
salt concentration. This prediction could be tested by computer simulations with explicit 
co- and counterions,— ^ and could be used to study experimentally the validity of the DLVO 
theory for oppositely charged colloids. 

Finally, we have studied the gas-liquid interfacial tension using histogram reweighting 
methods. We find that the dimensionless tension decreases for decreasing screening parame- 
ter. Upon scaling the interfacial tension with the contact value of the Coulomb interaction, 
we observed a collapse of the interfacial tensions onto a single curve. This means that for 
state points at coexistence and at the same scaled temperature Tq/T^ ct , the interfacial 
tension is determined solely by the contact value and not by the range of the interaction. 
There might be a possible connection with the well known similarities between the structures 
of the RPM and of the YRPM. Larsen and Rogde^ noted that their Monte Carlo results 
for the radial distribution functions of the YRPM were similar to those obtained for the 
RPM at different state points. The structurally equivalent states were further investigated 
by Copestake and Evans^ and by Leote de Carvalho and Evans 20 . This correspondence is 
due to the screening of charges in the RPM, which in turn is due to charge ordering . 20 ^^^ 
Consequently, the potential of mean force between two ions decays more rapidly than the 
bare Coulomb potential. For certain state points the potentials of mean force will be similar 
for the YRPM and the RPM. Preliminary results suggest that a similar explanation holds 
for the state points with similar interfacial tensions. Further investigation is in progress. 
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